home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
CD ROM Paradise Collection 4
/
CD ROM Paradise Collection 4 1995 Nov.iso
/
science
/
earthqua.zip
/
EQSIM
/
EQSIM.CPP
next >
Wrap
C/C++ Source or Header
|
1995-02-04
|
20KB
|
901 lines
// EQSIM.CPP - The simulator for the earthquake damage prevention program.
// Generates the building movement and displays it.
// Created by Misha Koshelev.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <dos.h>
#include <string.h>
#include <alloc.h>
#include <graphics.h>
#include "svga256.h"
#include "gendisp.h"
// Misc. definitions for init_globals
#define WAVELET_NUM 15
#define DELTA_OMEGA 0.2
#define FRACTAL_DIM 0.35
// Definitions for active and hybrid control
#define DELAY 3
#define CONTROL_NORMAL 2.0
#define HYBRID_SPOINT 1.0
// Program defaults
#define FELASTIC 1
#define FFRICTION 1
#define TSTEP 0.1
#define SNUM 300
#define BMASS 10.00
#define MAGNIFICATION 1.0
// Some defines
#define EXIT 255.255
typedef double (* cntrlfunc)(double); // The definition for the
// control function.
typedef struct
{
double f, a, v, x, e; // The force, accelaration,
// veloctity, energy used,
// and x movement of a
// building
// where the original position
// is 0.0
} bms;
bms far *blmov = NULL; // The building movement
// array.
double felastic_cos, ffriction_cos; // The coefficients for
// elastic and friction
// forces.
double tstep; // The time step used for
// the building movement.
int snum; // The number of steps.
double bmass; // The mass of the building.
// Predefined headers
int huge Svga256Detect(void);
double no_cntrl(double);
double active_cntrl(double);
double hybrid_cntrl(double);
double execute_cntrl(cntrlfunc, double);
void init_qglobals(double, double, double, double, int, double);
void init_gen_qglobals(double, double, double, int, double);
void done_qglobals(void);
void done_gen_qglobals(void);
void gen_blmov(cntrlfunc);
void plot_blmov(int, int, int, int, int, int, int);
void plot_disp(int, int, int, int, int, int, int);
void loadpal(char *);
void hrainbow(int, int, int, int, int, int, int, int);
void vrainbow(int, int, int, int, int, int, int, int);
double total_blmov(void);
double max_blmov(void);
double total_eused(void);
double show_control(cntrlfunc, char *, double);
void show_control_txt(cntrlfunc, char *, double);
double select_intensity(void);
// Global variables
int bpoly[10];
// Fill pattern, to fill the building with
char bpattern[8] = {0xFF, 0xFF, 195, 195, 195, 195, 0xFF, 0xFF};
// Set if program should write the result to a file (default)
int wtof = TRUE;
FILE far *tmf, far *mmf, far *ef;
// Set if demo mode (default off)
int demomode = FALSE;
// Set if no graphics (default off)
int graphics = TRUE;
// The main program.
//
#ifdef STANDALONE
void main(int argc, char *argv[])
#else
void eqsimmain(void)
#endif
{
double ti, i, intensity;
int ai;
int argprocessed = FALSE;
int gd = DETECT, gm, errorcode;
#ifdef STANDALONE
// Check for parameters
if (argc > 0)
{
for (ai = 1; ai < argc; ai++)
{
// If parameter is disable write to file, disable it
if (!strcmp(argv[ai], "/NOOUTPUTFILES"))
{
argprocessed = TRUE;
wtof = FALSE;
}
// If parameter is demo mode, enable it
if (!strcmp(argv[ai], "/DEMO"))
{
argprocessed = TRUE;
demomode = TRUE;
}
// If parameter is no graphics, disable it
if (!strcmp(argv[ai], "/NOGRAPHICS"))
{
argprocessed = TRUE;
wtof = TRUE;
demomode = FALSE;
graphics = FALSE;
}
if (argprocessed == FALSE)
{
printf("Error: Invalid parameter %s", argv[ai]);
exit(1);
}
argprocessed = FALSE;
}
}
#endif
#ifdef STANDALONE
if (graphics)
{
// Install the Super VGA 256 color driver
installuserdriver("Svga256", Svga256Detect);
// read the result of installation
errorcode = graphresult();
if (errorcode != grOk) // an error occurred
{
printf("Graphics error: %s\n", grapherrormsg(errorcode));
printf("Press any key to quit...");
getch();
exit(1);
}
// Initalize graphics mode
initgraph(&gd, &gm, "");
// read result of initialization
errorcode = graphresult();
if (errorcode != grOk) // an error occurred
{
printf("Graphics error: %s\n", grapherrormsg(errorcode));
printf("Press any key to quit...");
getch();
exit(1);
}
// Load the pallette
loadpal("std.pal");
}
#else
if (graphics)
cleardevice();
#endif
// If we should write to files, open the files
if (wtof)
{
tmf = fopen("totalmov.eqs", "wt");
mmf = fopen("maxmov.eqs", "wt");
ef = fopen("energy.eqs", "wt");
}
if (wtof && (!tmf || !mmf || !ef))
{
printf("Error: Could not create data files");
exit(1);
}
// Set up the TEX header
if (wtof)
{
fprintf(tmf,"\\settabs 4 \\columns\n\n");
fprintf(tmf,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
fprintf(mmf,"\\settabs 4 \\columns\n\n");
fprintf(mmf,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
fprintf(ef, "\\settabs 4 \\columns\n\n");
fprintf(ef,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
}
if (graphics)
{
// Initialize the polygon's coordinates
bpoly[0] = 300;
bpoly[1] = 249;
bpoly[2] = 340;
bpoly[3] = 249;
bpoly[5] = 200;
bpoly[7] = 200;
bpoly[8] = bpoly[0];
bpoly[9] = bpoly[1];
}
for (i=TSTEP; i<9.5; i+=TSTEP)
{
NewIntens:
// Write the intensity to the data files:
if (wtof)
{
fprintf(tmf, "\\+%.1f&", i);
fprintf(mmf, "\\+%.1f&", i);
fprintf(ef, "\\+%.1f&", i);
}
// No Control
intensity = pow(10, i/3);
init_qglobals(intensity, FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
if (!graphics)
{
show_control_txt(no_cntrl, "No Control:",i);
}
if (graphics)
{
ti = show_control(no_cntrl, "No Control:",i);
if (ti == EXIT)
goto Done;
if (ti != 0.0)
{
done_qglobals();
i = ti;
loadpal("std.pal");
goto NewIntens;
}
}
done_gen_qglobals();
// Show Active Control
init_gen_qglobals(FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
if (!graphics)
{
show_control_txt(active_cntrl, "Active Control:",i);
}
if (graphics)
{
ti = show_control(active_cntrl, "Active Control:",i);
if (ti == EXIT)
goto Done;
if (ti != 0.0)
{
done_qglobals();
i = ti;
loadpal("std.pal");
goto NewIntens;
}
}
done_gen_qglobals();
// Show Hybrid Control
init_gen_qglobals(FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
if (!graphics)
{
show_control_txt(hybrid_cntrl, "Hybrid Control:",i);
}
if (graphics)
{
ti = show_control(hybrid_cntrl, "Hybrid Control:",i);
if (ti == EXIT)
goto Done;
if (ti != 0.0)
{
done_qglobals();
i = ti;
loadpal("std.pal");
goto NewIntens;
}
}
if (wtof)
{
fprintf(tmf, "\\cr\n");
fprintf(mmf, "\\cr\n");
fprintf(ef, "\\cr\n");
}
done_qglobals();
}
if (wtof)
{
fclose(tmf);
fclose(mmf);
fclose(ef);
}
Done:
#ifdef STANDALONE
closegraph();
#endif
}
#ifdef STANDALONE
// Returns what video mode to use.
//
int huge Svga256Detect(void)
{
return SVGA640x480x256;
}
#endif
// The function for no control
//
// Input: Time
//
// Output: 0
//
double no_cntrl(double t)
{
return 0.0;
}
// The active control function.
//
// Input: Time
//
// Output: The number to add to the force
//
// The algorithm used is:
//
// f = -hx(t - T)
// control
//
// where h is a coefficient, t is the current time, and T depends on the
// delay of the control system.
//
double active_cntrl(double t)
{
double r = (-CONTROL_NORMAL * blmov[(int)(t / tstep) - DELAY].x);
blmov[(int)(t / tstep)].e = r * (blmov[(int)(t / tstep)].x - blmov[(int)
(t / tstep) - 1].x);
return r;
}
// The hybrid control function.
//
// Input: Time
//
// Output: The number to add to the force
//
// The algorithm used is:
//
// if x(t - T) <= p then:
//
// f = 0.0
// control
//
// else
//
// f = -hx(t)
// control
//
// where h is a coefficient, t is the current time,
// and T depends on the delay of the control system.
//
double hybrid_cntrl(double t)
{
double r = 0.0;
blmov[(int)(t/tstep)].e = 0.0;
if (blmov[(int)(t / tstep) - DELAY].x > HYBRID_SPOINT)
r = (-CONTROL_NORMAL * blmov[(int)(t / tstep)].x);
return r;
}
// Executes the control functions if the array element is greater than
// DELAY-1.
//
double execute_cntrl(cntrlfunc cf, double param)
{
if ((int)(param / tstep) - DELAY < 0)
{
blmov[(int)(param/tstep)].e = 0.0;
return 0.0;
}
return cf(param);
}
// Initializes the earthquake parameters.
//
// Input: Intensity, Elastic force coefficient, Friction force coefficient,
// the time step, the number of steps, and the building's mass.
//
// Output: None
//
void init_qglobals(double intens, double fe_cos, double ff_cos, double t_step,
int s_num, double b_mass)
{
snum = s_num;
if (blmov)
done_qglobals();
blmov = (bms far *)farcalloc(snum, sizeof(bms));
if (!blmov)
{
closegraph();
printf("Not enough memory\n");
exit(1);
}
for (int j=0; j<snum; j++)
{
blmov[j].x = 0.0; blmov[j].v = 0.0; blmov[j].f = 0.0; blmov[j].a = 0.0;
blmov[j].e = 0.0;
}
tstep = t_step;
felastic_cos = fe_cos;
ffriction_cos = ff_cos;
bmass = b_mass;
init_globals(intens, WAVELET_NUM, tstep * snum, DELTA_OMEGA, FRACTAL_DIM);
}
// Initializes the building movement parameters without initializing the
// earthquake parameters.
//
// Input: Elastic force coefficient, Friction force coefficient,
// the time step, the number of steps, and the building's mass.
//
// Output: None
//
void init_gen_qglobals(double fe_cos, double ff_cos, double t_step,
int s_num, double b_mass)
{
snum = s_num;
if (blmov)
done_gen_qglobals();
blmov = (bms far *)farcalloc(snum, sizeof(bms));
if (!blmov)
{
closegraph();
printf("Not enough memory\n");
exit(1);
}
for (int j=0; j<snum; j++)
{
blmov[j].x = 0.0; blmov[j].v = 0.0; blmov[j].f = 0.0; blmov[j].a = 0.0;
blmov[j].e = 0.0;
}
tstep = t_step;
felastic_cos = fe_cos;
ffriction_cos = ff_cos;
bmass = b_mass;
}
// Destroys the global variables
//
void done_qglobals(void)
{
farfree(blmov);
blmov = NULL;
done_globals();
}
// Destroys the global variables without destroying the earthquake parameters.
//
void done_gen_qglobals(void)
{
farfree(blmov);
blmov = NULL;
}
// Generates the building movement
//
// Input: A control function
//
// Output: None
//
// The following algorithm is used for every step where t is the current time.
//
// x(t) = v(t - s) * s + x(t - s)
//
// v(t) = a(t - s) * s + v(t - s)
//
// f(t) = f + f + f + f
// external elastic friction control
//
// f = -kx(t)
// elastic
//
// f = -nv(t)
// friction
//
// a(t) = f(t)
// ----
// m
//
// e(t) = f * x(t) - x(t - s)
// control
// /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\
// for active control and 0 for all others.
//
// where s is the time step, m is the mass of the building, k and n are
// coefficients, f(t) is the force in a current time, x(t) is the position
// related to the original position at the current time, v(t) is the velocity
// at the current time, and a(t) is the accelaration at the current time.
//
void gen_blmov(cntrlfunc cf)
{
// Generate values for time 0
blmov[0].x = 0.0;
blmov[0].v = 0.0;
blmov[0].f = disp(0.0) + execute_cntrl(cf, 0.0);
blmov[0].a = blmov[0].f / bmass;
// Generate all of the other values
for (int j=1; j<snum; j++)
{
blmov[j].x = blmov[j-1].v * tstep + blmov[j-1].x;
blmov[j].v = blmov[j-1].a * tstep + blmov[j-1].v;
blmov[j].f = disp(j * tstep) + (-felastic_cos * blmov[j-1].x) +
(-ffriction_cos * blmov[j-1].v) +
execute_cntrl(cf, j * tstep);
blmov[j].a = blmov[j].f / bmass;
}
}
// Plots the building movement.
//
void plot_blmov(int x, int y, int x1, int y1, int borcolor, int normcolor,
int limitcolor)
{
int cx, cy, middle, color;
int limit;
setcolor(borcolor);
rectangle(x, y, x1, y1); // Draw border rectangle
middle = ((y1-y) / 2) + y;
limit = x1+1;
if ((x1+1)-x > snum)
limit = x+snum;
for (cx=x;cx<limit;cx++)
{
cy = (int)blmov[cx-x].x * MAGNIFICATION + middle;
color = normcolor;
if (cy < y)
{
cy = y;
color = limitcolor;
}
if (cy > y1)
{
cy = y1;
color = limitcolor;
}
putpixel(cx, cy, color);
}
}
// Plots the earthquake displacements
//
void plot_disp(int x, int y, int x1, int y1, int borcolor, int normcolor,
int limitcolor)
{
int cx, cy, middle, color;
int limit;
setcolor(borcolor);
rectangle(x, y, x1, y1); // Draw border rectangle
middle = ((y1-y) / 2) + y;
limit = x1+1;
if (x1-x > snum)
limit = x+snum;
for (cx=x;cx<limit;cx++)
{
cy = (int)disp((cx - x) * tstep) + middle;
color = normcolor;
if (cy < y)
{
cy = y;
color = limitcolor;
}
if (cy > y1)
{
cy = y1;
color = limitcolor;
}
putpixel(cx, cy, color);
}
}
#ifdef STANDALONE
// Loads the palette from a raw palette file.
//
void loadpal(char *fn)
{
int color;
int r, g, b;
FILE *pal;
pal = fopen(fn, "rb");
if (pal == NULL)
{
closegraph();
printf("Error: Could not open palette file %s", fn);
exit(1);
}
for (color = 0; color < 256; color++)
{
fread(&r, 1, 1, pal);
fread(&g, 1, 1, pal);
fread(&b, 1, 1, pal);
setrgbpalette(color, r, g, b);
}
fclose(pal);
};
#endif
// Calculates the total building movements
//
double total_blmov(void)
{
double sum = 0.0;
for (int j=0; j<snum; j++)
sum += fabs(blmov[j].x);
return sum;
}
// Returns the maximum building movement.
//
double max_blmov(void)
{
double max;
if (snum == 1)
return fabs((double)blmov[0].x);
max = (fabs(blmov[0].x) > fabs(blmov[1].x)) ?
fabs(blmov[0].x)
: fabs(blmov[1].x);
if (snum == 2)
return max;
for (int j=2; j<snum; j++)
max = (max > fabs(blmov[j].x)) ? max : fabs(blmov[j].x);
return max;
}
// Returns the total energy used.
//
double total_eused(void)
{
double sum = 0.0;
for (int j=0; j<snum; j++)
sum += fabs(blmov[j].e);
return sum;
}
// Shows the control wave and the earthquake wave.
// Returns 0.0 or to what number to change the current intensity.
//
double show_control(cntrlfunc cf, char *cfname, double intens)
{
char tmp[20];
char ch;
int y;
// Set the colors
setfillstyle(SOLID_FILL, BLUE);
bar(1,1,640,300);
setfillstyle(SOLID_FILL, 79);
bar(1,250,640,300);
// If in demo mode, write the word DEMO
if (demomode)
{
sprintf(tmp, "DEMO");
settextstyle(TRIPLEX_FONT, HORIZ_DIR, 8);
setcolor(32);
outtextxy(300, 10, tmp);
setcolor(WHITE);
settextstyle(DEFAULT_FONT, HORIZ_DIR, 1);
}
// Show the intensity
sprintf(tmp, "Intensity: %.1lf", intens);
outtextxy(10, 10, tmp);
y = 10 + textheight(tmp) + 2;
// Show the earthquake displacement
sprintf(tmp, "Earthquake Disp:");
outtextxy(130, 300, tmp);
plot_disp(10,310,310,479,WHITE,WHITE,RED);
// Show the building movement
gen_blmov(cf);
outtextxy(440, 300, cfname);
plot_blmov(320,310,620,479,WHITE,WHITE,RED);
// Display the data
sprintf(tmp, "Total: %.3lf", total_blmov());
outtextxy(10,y,tmp);
if (wtof)
fprintf(tmf, "%.3lf&", total_blmov());
y += textheight(tmp);
sprintf(tmp, "Max. Move.: %.3lf", max_blmov());
outtextxy(10,y,tmp);
if (wtof)
fprintf(mmf, "%.3lf&", max_blmov());
y += textheight(tmp);
sprintf(tmp, "E. Used: %.3lf", total_eused());
outtextxy(10,y,tmp);
if (wtof)
fprintf(ef, "%.3lf&", total_eused());
y += textheight(tmp);
bpoly[4] = 340;
bpoly[6] = 300;
setfillstyle(SOLID_FILL, 128);
fillpoly(5, bpoly);
setfillpattern(bpattern, WHITE);
fillpoly(5, bpoly);
setfillstyle(SOLID_FILL, BLUE);
for (int k=0; k<SNUM; k++)
{
if (bpoly[4] != (int)blmov[k].x + 340 && bpoly[6] !=
(int)blmov[k].x + 300)
{
bar(1,y+1,640, 249);
bpoly[4] = (int)blmov[k].x + 340;
bpoly[6] = (int)blmov[k].x + 300;
setfillstyle(SOLID_FILL, 128);
fillpoly(5, bpoly);
setfillpattern(bpattern, WHITE);
fillpoly(5, bpoly);
setfillstyle(SOLID_FILL, 124);
setfillstyle(SOLID_FILL, BLUE);
}
else
delay(17); // If no time taken to draw delay for 17 ms
}
if (!demomode)
{
#ifdef STANDALONE
outtextxy(200, 265, "Press any key to continue...");
outtextxy(200, 275, "To change intensities press <I>");
outtextxy(200, 285, "To exit press <ESC>");
#else
outtextxy(200, 265, "Press any key to continue...");
outtextxy(200, 275, "To exit press <ESC>");
#endif
while (!kbhit());
ch = getch();
setfillstyle(SOLID_FILL, BLACK);
bar(1,1,640,480);
#ifdef STANDALONE
if (ch == 'i' || ch == 'I')
return select_intensity();
#endif
if (ch == 27)
{
done_qglobals();
if (wtof)
{
fclose(tmf);
fclose(mmf);
fclose(ef);
}
#ifdef STANDALONE
closegraph();
#endif
return EXIT;
}
}
if (demomode)
{
if (kbhit())
{
done_qglobals();
if (wtof)
{
fclose(tmf);
fclose(mmf);
fclose(ef);
}
#ifdef STANDALONE
closegraph();
#endif
return EXIT;
}
setfillstyle(SOLID_FILL, BLACK);
bar(1,1,640,480);
}
return 0.0;
}
// Shows the current intensity, etc. and writes the data to a file.
//
void show_control_txt(cntrlfunc cf, char *cfname, double intens)
{
char tmp[20];
char ch;
// Show the intensity
sprintf(tmp, "Intensity: %.1lf", intens);
printf("%s", tmp);
printf("\n");
printf("%s\n", cfname);
cprintf("Please wait...");
// Generate the building movement
gen_blmov(cf);
// Display the data
sprintf(tmp, "Total: %.3lf", total_blmov());
printf("\n%s\n", tmp);
if (wtof)
fprintf(tmf, "%.3lf&", total_blmov());
sprintf(tmp, "Max. Move.: %.3lf", max_blmov());
printf("%s\n", tmp);
if (wtof)
fprintf(mmf, "%.3lf&", max_blmov());
sprintf(tmp, "E. Used: %.3lf", total_eused());
printf("%s\n\n", tmp);
if (wtof)
fprintf(ef, "%.3lf&", total_eused());
if (kbhit())
{
ch = getch();
done_qglobals();
if (wtof)
{
fclose(tmf);
fclose(mmf);
fclose(ef);
}
exit(0);
}
}
#ifdef STANDALONE
// Asks the user for a new intensity.
//
double select_intensity(void)
{
double r;
int gd = DETECT, gm;
restorecrtmode();
printf("What intensity do you want to switch to (X.X)? ");
scanf("%3lf", &r);
if (r > 9.5)
r = 9.5;
setgraphmode(getgraphmode());
return (r);
}
#endif